#include <cilk/cilk.h>
#include <iostream>
#include <iomanip>
#include <limits>
#include <time.h> 
#include <string>
#include <fstream>

#include "../rungekutta2_serial/rungekutta2.h"
#include "../rungekutta2_serial/gtest.h"
#include "../rungekutta2_serial/Timer.h"

using namespace std;

void cilkplus_midpointMethod(double *a, size_t n, double h, double *y0, double *res, int steps);


int main(int argc, char **argv)
{
	double *pA, h, *y0, *res;
	ofstream out("log.txt");
	Stopwatch *s = createStopwatch();
	for (size_t n = 2, steps = n * 128; n <= 1<<9; n*=2, steps = n * 128)
	{
		s->reset();
		
		pA = new double[n * n];
		y0 = new double[n];
		res = new double[n * steps];
		RandomInitData(pA, h, y0, n);
		//alg begin
		s->start();
		cilkplus_midpointMethod(pA, n, h, y0, res, steps);
		s->stop();
		//alg end
		out << "{" << n << "," << s->getElapsed() << "},";
		if (!IsCorrectResult(pA, n, h, y0, res, steps))
			printf("error n=%d", n);
		delete[] pA;
		delete[] y0;
		delete[] res;
	}

	return 0;
}

double dotProd(const double *a, const double *b, const int n) 
{ 
	double res = 0.0; 
	for (int i = 0; i < n; i++) 
	{ 
		res += a[i] * b[i]; 
	}
	return res; 
} 

void cilkplus_MatrixVectorMult(double* pMatrix, double* pVector, double* pResult, int n, int m) 
{ 
	cilk_for (int i = 0; i < n; i++) 
	{ 
		pResult[i] = dotProd(&(pMatrix[i * m]), pVector, m); 
	} 
}

/*
*Midpoint method is method Runge-Kutta second-order:
*x(0)
*x(n+1) = x(n)+h
*y(n+1) = y(n) + h*f(  x(n) + h/2,    y(n) + h/2*f(x(n), y(n))  )
*
*Differential equation:
*u'=Au   A - Rnxn
*f(x,y) = Ay
*
*Algorithm:
*y(n+1) = y(n) + h*f(  x(n) + h/2,    y(n) + h/2*f(x(n), y(n))  ) =
*= y(n) + h*A(y(n) + h/2*A*y(n)) = 
*= (E + h*A + h*h/2*A*A)y(n) = 
*= Zy(n)
*
*y(n+1) = Zy(n)
*/
void cilkplus_midpointMethod(double *a, size_t n, double h, double *y0, double *res, int steps)
{
	__cilkrts_set_param("nworkers", "2");
	double *z = new double[n * n];
	matrixMatrixMult(a, a, z, n, n, n);
	matrixScalarMult(z, n, n, h/2);
	matrixAdd(z, a, n, n);
	matrixScalarMult(z, n, n, h);
	matrixEAdd(z, n);
	copy(y0, y0 + n, res);
	for(int i = 0; i < steps - 1; i++)
	{
		cilkplus_MatrixVectorMult(z, res + i * n, res + (i + 1)*n, n, n);
	}

	delete [] z;
}